/*! 
 ***************************************************************************
 * \file   convection3.cpp
 * Forced Convection-Diffusion in 3D.
 *
 ***************************************************************************
 * \author    Luis M. de la Cruz [ Fri Jun  8 12:03:41 CDT 2007 ]
 ***************************************************************************
 * This code solves the next equation:
 * \f[ 
 * \frac{\partial T}{\partial t} + 
 * u_j \frac{\partial T}{\partial x_j} = \nabla^2 T
 * \,\,\,\,\,\, for \,\,\,\,\, j=1,2,3
 * \f]. 
 * \par
 * This equation is solved in a cubic prism \f$ x, y, z \in [0,1] \f$. 
 * The velocity is prescribed and has the form:
 * \f[
 * u_1 = -A \cos(\pi y) \sin(\pi \lambda x / L_x );
 * \f]
 * \f[
 * u_2 = \frac{A \lambda}{L_x} \sin(\pi y) \cos(\pi \lambda x / L_x);
 * \f]
 * \f[
 * u_3 = 0
 * \f]
 * The boundary conditions are as shown in the next figure:
 *
 * \image html convection3.png "Domain in 3D" width=10cm
 * \image latex convection3.eps "Domain in 3D" width=10cm
 *
 * The equation can be solved using Upwind, Central Average (CDS) or QUICK 
 * schemes for the approximation of convective terms. CDS and QUICK can
 * be used with \c CDS_CoDi and \c QUICK_CoDi, or with \c CDS_Hay and 
 * \c QUICK_Hay. The latter two are the deferred implementation from
 * Hayase et al. 1992, which are more stable.
 *
 * The results are as presented in the next figure:
 *
 * \image html  convForced3D01.png "Final result" width=10cm 
 * \image latex  convForced3D01.eps "Final result" width=10cm 
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 * \par
   \verbatim
   % make
   % ./convection3 \endverbatim
 * 
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"

#ifdef WITH_CDS_CoDi
#include "Schemes/CDS_CoDi.hpp"
#endif

#ifdef WITH_CDS_Hay
#include "Schemes/CDS_Hay.hpp"
#endif

#ifdef WITH_Upwind_CoDi
#include "Schemes/Upwind_CoDi.hpp"
#endif

#ifdef WITH_Quick_CoDi
#include "Schemes/Quick_CoDi.hpp"
#endif

#ifdef WITH_Quick_Hay
#include "Schemes/Quick_Hay.hpp"
#endif

#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,3>::huge ScalarField3D;

template<class sfield>
void initialVelocity(sfield &u, sfield &v, sfield &w, double dx, double dy,
		      double A, double l_x, double l)
{
    const double PI = 3.14159265359;
    const int bi = u.lbound(firstDim) + 1, ei = u.ubound(firstDim) - 1,
	bj = u.lbound(secondDim) + 1, ej = u.ubound(secondDim) - 1,
	bk = u.lbound(thirdDim) + 1, ek = u.ubound(thirdDim) - 1;

    double A2 = A * l / l_x;

    for(int k = bk; k <= ek; ++k) 
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j) {
		u(i,j,k) = -A * cos (PI * dy * (j-0.0)) * 
		    sin (PI * l * dx * (i-0.0) / l_x );    
		v(i,j,k) = A2 * sin (PI * dy * (j-0.0)) * 
		    cos (PI * l * dx * (i-0.0) / l_x);
		w(i,j,k) = 0.0;
	    }
}

int main() 
{
    double length_x, length_y, length_z;    
    int num_nodes_x, num_nodes_y, num_nodes_z;
    double left_wall, right_wall;
    double tolerance, dt, error = 1.0, residual;
    double lambda, amplitude;
    int max_iter, nmax, frequency, iteration = 0;

// Reading data from the input file
    std::ifstream input_cin ("input3");
    input_cin >> length_x
	      >> length_y
	      >> length_z
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> num_nodes_z
	      >> right_wall 
	      >> left_wall 
	      >> dt
	      >> nmax
	      >> tolerance
	      >> max_iter
	      >> lambda
	      >> amplitude
	      >> frequency;
    input_cin.close();

    StructuredMesh<Uniform<double, 3> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y,
					     length_z, num_nodes_z);
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    double dz = mesh.getDelta(Z);
    mesh.print();
    int num_vols_x = mesh.getExtentVolumes(X);
    int num_vols_y = mesh.getExtentVolumes(Y);
    int num_vols_z = mesh.getExtentVolumes(Z);

    ScalarField3D T(mesh.getExtentVolumes());
    ScalarField3D u(mesh.getExtentVolumes());
    ScalarField3D v(mesh.getExtentVolumes());
    ScalarField3D w(mesh.getExtentVolumes());
    ScalarField3D us(num_nodes_x, num_vols_y, num_vols_z); // staggered u-velocity 
    ScalarField3D vs(num_vols_x, num_nodes_y, num_vols_z); // staggered v-velocity 
    ScalarField3D ws(num_vols_x, num_vols_y, num_nodes_z); // staggered w-velocity 

    ScalarField3D Tn(mesh.getExtentNodes());
    ScalarField3D un(mesh.getExtentNodes());
    ScalarField3D vn(mesh.getExtentNodes());
    ScalarField3D wn(mesh.getExtentNodes());

    Range all = Range::all();
    T(T.lbound(firstDim), all, all) = left_wall;   
    T(T.ubound(firstDim), all, all) = right_wall;  

    initialVelocity(un, vn, wn, dx, dy, amplitude, length_x, lambda);
    NumUtils::staggerFromNodesX(us, un);
    NumUtils::staggerFromNodesY(vs, vn);
    NumUtils::staggerFromNodesZ(ws, wn);

#ifdef WITH_DX
    InOut::writeToFile_DX(un, vn, wn, 0, "./Data3/velc.", dx, dy, dz);
#endif

    SparseMatrix< Diagonal<double, 3 > > A(num_nodes_x, num_nodes_y, num_nodes_z); 

    ScalarField3D              b(num_nodes_x, num_nodes_y, num_nodes_z);

#ifdef WITH_CDS_CoDi
    ScalarEquation<CDS_CoDi<double,3> > conv_diff(T, A, b, mesh.getDeltas()); 
#endif

#ifdef WITH_CDS_Hay
    ScalarEquation<CDS_Hay<double,3> > conv_diff(T, A, b, mesh.getDeltas()); 
#endif

#ifdef WITH_Upwind_CoDi
    ScalarEquation<Upwind_CoDi<double,3> > conv_diff(T, A, b, mesh.getDeltas());
#endif

#ifdef WITH_Quick_CoDi
    ScalarEquation<Quick_CoDi<double,3> > conv_diff(T, A, b, mesh.getDeltas());    
#endif

#ifdef WITH_Quick_Hay
    ScalarEquation<Quick_Hay<double,3> > conv_diff(T, A, b, mesh.getDeltas());    
#endif

    conv_diff.setDeltaTime(dt);
    conv_diff.setGamma(1.0);
    conv_diff.setNeumann(TOP_WALL);
    conv_diff.setNeumann(BOTTOM_WALL);
    conv_diff.setDirichlet(LEFT_WALL);
    conv_diff.setDirichlet(RIGHT_WALL);
    conv_diff.setNeumann(FRONT_WALL);
    conv_diff.setNeumann(BACK_WALL);
    conv_diff.setUvelocity(us);
    conv_diff.setVvelocity(vs);
    conv_diff.setWvelocity(ws);
    conv_diff.print();

    ofstream error_fp("./Data3/error"), residual_fp("./Data3/residual"), 
      nusselt_fp("./Data3/nusselt");
    
    while ( (error > tolerance) && (iteration <= nmax) ) {
      conv_diff.calcCoefficients();
      Solver::TDMA3D(conv_diff, tolerance, max_iter);
      residual = conv_diff.calcResidual();
      error = NumUtils::calcErrorL2(T, conv_diff.phi);
      
      conv_diff.update();
      
      error_fp << iteration << "\t" << error << endl;
      //	residual_fp << iteration << "\t" << Solver::residual << endl;
      //	nusselt_fp << iteration  << "\t" 
      //		   << NumUtils::calcNusselt(T,dx,dy,dz,0) << endl;
      
      cout << "\n ---- Iteration = " << ++iteration 
	   << " Residual = " << residual
	   << " Error = " << error
	   << endl;
      
      if (!(iteration % frequency)) {
	NumUtils::interpolateToNodes(Tn, T);

#ifdef WITH_DX
	InOut::writeToFile_DX(Tn, iteration, "./Data3/temp.", dx, dy, dz);
#endif

#ifdef WITH_AM
	InOut::writeToFile_AM(Tn, iteration, "./Data3/temp.", dx, dy, dz);
#endif
      }
    }
    
#ifdef WITH_AM
    std::cout << "\n\n NOTE: The data were written in the [./Data2] directory in AmiraMesh format";
#endif

#ifdef WITH_DX
    std::cout << "\n\n NOTE: The data were written in the [./Data3] directory in DX format";
#endif

    return 0;
}







